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Abstract 

The spectrum and coherency are useful quantities for characterizing the 
temporal correlations and functional relations within and between point 
processes. This paper begins with a review of these quantities, their inter- 
pretation and how they may be estimated. A discussion of how to assess 
the statistical significance of features in these measures is included. In ad- 
dition, new work is presented which builds on the framework established 
in the review section. This work investigates how the estimates and their 
error bars are modified by finite sample sizes. Finite sample corrections 
are derived based on a doubly stochastic inhomogeneous Poisson process 
model in which the rate functions are drawn from a low variance Gaussian 
process. It is found that, in contrast to continuous processes, the variance 
of the estimators cannot be reduced by smoothing beyond a scale which is 
set by the number of point events in the interval. Alternatively, the degrees 
of freedom of the estimators can be thought of as bounded from above by 
the expected number of point events in the interval. Further new work de- 
scribing and illustrating a method for detecting the presence of a line in a 
point process spectrum is also presented, corresponding to the detection of 
a periodic modulation of the underlying rate. This work demonstrates that 
a known statistical test, applicable to continuous processes, applies, with 
little modification, to point process spectra, and is of utility in studying a 
point process driven by a continuous stimulus. While the material discussed 
is of general applicability to point processes attention will be confined to 
sequences of neuronal action potentials (spike trains) which were the moti- 
vation for this work. 
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1 Introduction 

The study of spike trains is of central importance to electrophysiology. Often 
changes in the mean firing rate are studied but there is increasing interest in 
characterising the temporal structure of spike trains, and the relationships between 
spike trains, more completely ( Gray et al., 198H ; Gerstein et al., 1985| ; Abclcs] 



|et al., 1983|) . A natural extension to estimating the rate of neuronal firing is to 
estimate the autocorrelation and the cross-correlation functions^. This paper will 
discuss the frequency domain counterparts of these quantities. Auto- and cross- 
correlations correspond to spectra and cross spectra respectively. The coherency, 
which is the normalised cross spectrum, does not in general have a simple time 
domain counterpart. 

The frequency domain has several advantages over the time domain. Firstly 
often subtle structure can be detected with the frequency domain estimators which 
is difficult to observe with the time domain estimators. Secondly, the time domain 
quantities have problems which are associated with sensitivity of the estimators 
to weak non-stationarity and the non-local nature of the error bars flBrody, 19981 ) • 
These problems are greatly reduced in the frequency domain. Thirdly, reasonably 
accurate confidence intervals may be placed on estimates of the second order 
properties in the frequency domain which permits the statistical significance of 
features to be assessed. Fourthly, the coherency provides a normalised measure of 
correlations between time series, in contrast with time domain cross-correlations 
which are not normalisable by any simple means. 

This paper begins by reviewing the population spectrum and coherency for 
point processes and motivating their use by describing some example applications. 
Next direct, lag window and multitaper estimators of the spectrum and coherency 
are presented. The concept of degrees of freedom is introduced and used to obtain 
large sample error bars for the estimators. Many elements of the work discussed 
in the review section of this paper can be found in the references ( [Percival and 



Walden, 1993| |Cox and Lewis, 19661 ; [Brillinger, 19781 ; [Bartlett, 19661) . Most of the 



material in these references is targeted at either spectral analysis of continuous 
processes or at the analysis of point processes but with less emphasis on spectral 
analysis. Building on this framework corrections, based on a specific model, will 
be given for finite sample sizes. These corrections are cast in terms of a reduction 
in the degrees of freedom of the estimators. For a homogeneous Poisson process 
the modified degrees of freedom is the harmonic sum of the the asymptotic de- 



1 Definitions of these quantities will be given in section 2.3 
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grees of freedom and twice the number of spikes used to construct the estimate. 
Modifications to this basic result are given for structured spectra and tapered 
data. A section is included on the treatment of point process spectra which con- 
tain lines. A statistical test for the presence of a line in a background of coloured 
noise is given, and the method for removal of such a line described. An example 
application to periodic stimulation is given. 

2 Population measures and their interpretation 



2.1 Counting representation of a spike train 

A spike train may be regarded as a point process. If the spike shapes are neglected, 
it is completely specified by a series of spike times and the start and end points 
of the recording interval [0, T\. It is convenient to introduce some notation which 
enables the subsequent formulae to be written in a compact form ( [Brillinger^ 
1978 ). The counting process N(t) is defined as the number of spikes which occur 
between the start of the interval (t = 0) and time t. The counting process has the 
property that the area beneath it grows as t becomes larger. This is undesirable 
because it leads to an interval dependent peak at low frequencies in the spectrum. 
To avoid this a process N(t) = N(t) — Xt, where A is the mean rate, which has zero 
mean may be constructed. Note that dN(t) = N(t + dt) — N(t) which is either 
1 — Xdt or — Xdt depending on whether or not there is a spike in the interval 
dt. Thus dN(t)/dt is a series of delta functions^ with the mean rate subtracted. 
Figure [I] illustrates the relationship between N(t), N(t), and dN(t)/dt. 



N(t) 
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Figure 1: Example illustrating how the processes N, N and dN /dt relate to each 
other. The vertical lines in the process dN /dt depict delta functions. 



2 A delta function is a generalized function. It has an area of one beneath it but has zero 
width and therefore infinite height. 
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2.2 Stationarity 



It will be assumed in what follows that the spike trains are second order stationary. 
This means that their first and second moments do not depend on the absolute 
time. In many electrophysiology experiments this is not the case. In awake 
behaving studies the animal is often trained to perform a highly structured task. 
Nevertheless it may still be the case that over an appropriately chosen short 
time window, the statistical properties are changing slowly enough for reasonable 
estimates of the spectrum and coherency to be obtained. As an example, neurons 
in primate parietal area PRR exhibit what is known as memory activity during 
a delayed reach task ( jSnyder et al., ICW7 ). The mean firing rate of these neurons 
varies considerably during the task but during the memory period is roughly 
constant. The assumption of stationarity during the memory period is equivalent 
to the intuitive notion that there is nothing special about 0.75s into the memory 
period as compared to say 0.5s. Second order stationarity implies that the mean 
firing rate (A) is independent of time and additionally that the autocovariance 
depends only on the lag (r) and not on the absolute time. 



2.3 Definitions 

Equations [I] - |] give the first and second 
for a stationary process. The spectrum 
autocovariance function (/x(r) + X8(t)). 



order moments of a single spike train 
S(f) is the Fourier transform of the 



E{dN(t)} _ 

df A (1) 

^P=0 (2) 

,J +xs(T)= Emm±iA (3) 

/CO 

fi(T)exp(-27rifr)dT (4) 
-oo 

Where E denotes the expectation operator. 



The autocovariance measures how likely it is that a spike will occur at time t+r 
given that one has occurred at time t. Usually /i(r) is estimated rather than the 
full autocovariance which includes a delta function at zero lagQ. However, in order 
to take the Fourier transform the full autocovariance is required. The inclusion 
of this delta function leads to a constant offset of the spectrum. This offset is an 

3 When estimating the autocovariance using a histogram method one generally omits the 
spike at the start of the interval which would always fall in the bin nearest zero. 
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important difference between continuous time processes and point processes. The 
population coherency j(f) is defined in equations |5| - [7| 



. E[dN a (t)dN b (t + r)} 

^ b{r) = m (5) 

/oo 
/jL ab (r) exp(-2vrz fr)dr + X a 5 ab (6) 
-oo 

1(f) = . Sl2{f) (7) 
pix{f)S 22 {f) 

Where indices 1 and 2 denote simultaneously recorded spike trains from different 
cells. 

Unlike the spectrum, which is strictly real and positive, the coherency is a 
complex quantity. The modulus of the coherency, which is known as the coher- 
encef], can only vary between zero and one. This makes coherence particularly 
attractive for detecting relationships between spike trains as it is insensitive to 
the mean spike rates. 



3 Examples and their interpretation 

Before discussing the details regarding how to estimate the spectrum and co- 
herency it will be helpful to motivate them further by considering some simple 
examples. 

3.1 Example population spectra 

For a homogeneous Poisson process of constant rate A the autocovariance is simply 
X5(t) and hence the spectrum is a constant equal to the rate A. At the opposite 
extreme consider the case where the spikes are spaced by intervals At. This is not 
a stationary process but if a small amount of drift is permitted, so that over an 
extended period there is nothing special about a given time, it becomes stationary. 
The spectrum of this process contains sharp lines at integer multiples of / = 
Due to the drift the higher harmonics will become increasingly blurred and in the 
high frequency limit the spectrum will tend towards a constant value of the mean 
rate A. As a final example consider the case where /x(r) is a negative Gaussian 
centered on zero r. This form of fi{r) is consistent with the probability of firing 
being suppressed after firing^. The spectrum of this process will be below A at 

4 Some authors define coherence as the modulus squared of the coherency. 
5 This need not necessarily correspond to the biophysical refractive period but, it could arise, 
rather from a characteristic integration time. 
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low frequencies and will go to a constant value A at high frequencies. Figure |2| 
illustrates these different population spectra. 



S(f) 



S(f) 
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X 




(a) f (b) f (c) f 

Figure 2: Example population spectra for different types of underlying process, 
(a) Homogeneous Poisson process with rate A. (b) Regularly spaced spikes with 
jitter, (c) Spike trains in which the probability of firing is suppressed immediately 
after firing. 



3.2 Example population coherency 

The population coherency of two homogeneous Poisson processes is zero. In con- 
trast if two spike trains are equal then the coherence is one and the phase of the 
coherency is zero at all frequencies. If two spike trains are identical but offset by 
a lag At then the coherence will again be one but the phase of the coherency 
will vary linearly with frequency with a slope proportional to At and given by 
o(f) i^/Av. 

4 Estimating the spectrum 

Section [3| demonstrated that the population spectrum may provide insights into 
the nature of a spike train. In this section the question of how to estimate the 
spectrum from a finite section of data will be introduced. In what follows the 
population quantity A in the definition of N(t) is replaced by a sample estimate 
N(T)/T. 

4.1 Direct Spectral Estimators 
4.1.1 Definition 

A popular, though seriously flawed, method for estimating the spectrum is to take 
the modulus squared of the Fourier transform of the data dN(t). This estimate 
is known as the Periodogram and is the simplest example of a direct spectral 
estimator. More generally, a direct spectral estimator is the modulus squared 
of the Fourier transform of the data but with the data being multiplied by an 
envelope function h(t), known as a taper ( [Percival and Walden, 1993|) . Equations 
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|H] define the direct estimator. On substituting N(t) into equation || a form 



amenable to implementation on a computer is obtained (equation [II]). In this 
form the Fourier transform may be computed rapidly and without the need for 
the binning of data. Note that equation [H^ results in h(t) scaling as l/y/T as 
the sample length is altered. This ensures proper normalization of the Fourier 
transformation as sample size varies. 

I D (f) = \J D (f)\ 2 (8) 
J D (f)= [ T h(t)e- 2mft dN(t) (9) 

Where, 



o 



o 



T 

hitfdt = 1 (10) 
J°(f) = £ h(t 3 )e-^ - iV(T ^ (/) (11) 



and H(f) is the Fourier transform of the taper. 



The direct estimator suffers from bias and variance problems, described below, 
and is of no practical relevance for a single spike train sample. 



4.1.2 Bias 

It may not be immediately apparent why the above procedure is an estimate of the 
spectrum, especially when one is permitted to multiply the data by an arbitrary, 
albeit normalized, taper. The relation between I D {f) and the spectrum may be 
obtained by taking the expectation of equation ||. 



/OO POO 
/ h^h^e-^^-'UN^dNit')} (12) 
-oo J —oo 

Assuming that the integration and expectation operations may be interchanged 
and substituting equation ^| yields 



/oo roo 
/ /i(t)/i(t>- 2m/ ('-' >{//(* - t') + X8(t - t')}dtdt' (13) 
-oo J ~oo 

Which may be rewritten in the Fourier domain as, 

/oo 
S{f)\HU - f)\ 2 df (14) 
-oo 

6 For the moment, we assume that the population quantity A is known. This is of course not 
the case in practice, and one employs the estimate N(T)/T as stated before. The effect of this 
extra uncertainty is given in equation Hq. 
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The expected value of the direct estimator is a convolution of the true spectrum 
and the modulus squared of the Fourier transform of the taper. The normaliza- 



tion condition on the taper (equation 1C) is equivalent to the requirement that the 
kernel of the convolution has unit area underneath it. Sharp features in the true 
spectrum will be thus be smeared by an amount which depends on the width of 
the taper in the frequency domain. If the taper is well localized in the frequency 
domain the expected value of the direct estimate is close to the true spectrum but 
if the taper is poorly localized then the expected value of the direct estimator will 
be incorrect i.e. the direct spectral estimator is biased. There is a fundamental 
level beyond which the bias cannot be reduced, due to the uncertainty relation for- 
bidding simultaneous localization of a function in the time and frequency domains 
below a given limit. Since the maximum width of the taper is T the minimum 
frequency spread is 1/T which is known as the Raleigh frequency. Figure [3] shows 
the smoothing kernel for a rectangular taper and a T of 0.5s. Note that this kernel 
has large sidelobes which is the primary motivation for using tapering. 




frequency (Hz) 

Figure 3: The smoothing kernel \H(f)\ 2 . This is the expected direct estimate 
of the spectrum in the case of a population spectrum which has a delta function 
(very sharp feature) at the center frequency. A rectangular taper of length 0.5s 
was used. Solid vertical lines are drawn at ± the Raleigh frequency. 

In the above argument equation || was used in spite of the appearance of the 
population quantity A rather than the sample estimate N(T)/T for which equation 
|12| was defined. A more careful treatment, which includes this correction, leads to 
an additional term at finite sample sizes in the expectation of the direct spectral 
estimator at low frequencies. The full expression is given below, 



/CO 
S(f)\H(f-f')\ 2 df-\H(f)\ 2 S(0)/T 
-co 



(15) 



In the case of the periodogram, where h{t) = 1/vT, the effect is clear since in 
this case J D (0) = and hence / D (0) = for any set of spike times and any T. 



4.1.3 Asymptotic variance 

In the previous section it was shown that provided the taper is sufficiently local in 
frequency the expected value of the direct spectral estimator will be close to the 



true spectrum. However, the fact that the estimate is on average close to the true 
spectrum belies a serious problem with direct spectral estimators, namely that the 
estimates have very large fluctuations about this mean. The underlying source of 
this problem is that one is attempting to estimate the value of a function at an 
infinite number of points using a finite sample of data. The problem manifests 
itself in the fact that direct spectral estimators are inconsistent estimators of the 
spectrumP]. In fact it may be shown that, under fairly general assumptions, the 
estimates are distributed exponentially (or equivalently as S(f)xz/2) for asymp- 
totic sample sizes (i.e. T — > oo) QBrillinger, 1972|) . Figure |] illustrates that 
direct spectral estimators are noisy and untrustworthy, a fact emphasised by the 
observation that the x\ distribution has a standard deviation equal to its mean. 
In the next three subsections methods for reducing the variance of direct spectral 
estimators using different forms of averaging will be discussed. 

moo 

100 

E 1" 

C 1 
a> 

Q. 

W 0.1 

0.01 



frequency (Hz) 

Figure 4: An example of a direct spectral estimate. A 40% cosine taper was used. 
A sample of duration 20s was drawn from a homogeneous Poisson process with 
a constant rate of 50 Hz. The population spectrum for this process is flat and is 
shown by the solid horizontal line. The direct spectral estimate is clearly noisy 
although on average the correct spectrum is obtained. 

4.2 Trial averaging 

If there are a number of trials (iV r ) available then the variance of the direct 
estimator may be reduced by trial averaging. 

J OT (/)= ArEtfC/) ( 16 ) 

n= i 

Where I^if) is the direct spectral estimate based on the n th trial. 

In the large T limit taking the average entails summing Nt independent sam- 
ples from a x\ distribution the result of which is distributed as x\n t ■ The reduc- 
tion in variance is inversely proportional to the number of trials corresponding to 
a reduction in standard deviation which is the familiar factor of 1/y/Nr- 

7 Inconsistent estimators have a finite variance even for an infinite length sample. 
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At first sight it appears one would be getting something for nothing by breaking 
a single section of data into segments and treating them as separate trials. This 
is, of course, not the case. The reason is that if the data is segmented into short 
length samples, there is loss of frequency resolution proportional to the inverse of 
segment length. Lag window and multitaper estimators use the information from 
these independent estimates without artificially segmenting the data. 



4.3 Lag Window Estimates 



A powerful property of the frequency domain is that, unless two frequencies are 
very close together, direct estimates of the spectrum of a stationary process at 
different frequencies are nearly uncorrelated. This property arises when the covari- 
ance between frequencies falls off rapidly. If the true spectrum varies slowly over 
the width of the covariance then the large sample covariance of a direct spectral 



estimator is given by equation 17 



cov{I D (f 1 ),I D (f 2 )} ~ E{I D (f)} 2 
Where / = (/i + / 2 )/2 and Af — fx — f 2 



h(t) 2 e~ 2mAft dt 



(17) 



For Af = 0, this expression reduces to the previously mentioned result that 
the variance of the estimator is equal to the square of the mean. For Af » 1/T, 
I JZo h(t) 2 e- 2wiAft dt\ 2 -> 0, since h(t) 2 is a smooth function with extent T. This 
implies that cov{I D (fi), I D (f 2 )} ~ for \f 1 — f 2 \ » 1/T. The approximate 
independence of nearby points means that, if the true spectrum varies slowly 
enough, then closely spaced points will provide several independent estimates of 
the same underlying spectrum. This is the motivation for the lag window estimator 
which is simply a smoothed version of the direct spectral estimator ( [Percival and 



Walden, 1993|) . The lag window estimator is defined in equations pl| and [19 



/oo 
K(f - f)I D (f)df (18) 
-oo 

Where, 

poo 

K(f)df = 1 (19) 



Averaging over trials may be included by using the trial averaged direct spec- 
tral estimate I DT (see equation [16]) in place of I D in the above expression. It is 



assumed that K(f) is a smoothing kernel with reasonable properties. 



4.3.1 Bias 

The additional smoothing of the lag window kernel modifies the bias properties 
of the estimator from those expressed in equation |T3|. The expected value of the 
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lag window estimator is given by, 



E{I LW (f)} = [°° K(f-f')\H(f'-f")\ 2 S(f")df'df"-^- f°° K(f-f')\H(f')\ 2 df 

J —oo 1 J —oo 

(20) 

4.3.2 Asymptotic Variance 



The large sample variance of this estimator is readily obtained using equation 17 



var{I LW (f)} = ^E{I LW (f)} 2 (21) 

/oo poo 
/ K(f)K(f)\H(f - f)\ 2 dfdf (22) 
-oo J —oo 



H(f) = / h{t) 2 e- 27rift dt (23) 



Where, 



and, 



Equation [21] includes the reduction in variance due to trial averaging. l/£ 
can be interpreted as the effective number of independent estimates beneath 
the smoothing kernel, as demonstrated by the following qualitative argument. 
If Af is the frequency width of the smoothing kernel K(f) and Sf is the fre- 
quency width of the taper TC(f) then since K(f) ~ 1/A/ it follows that £ ~ 
1/(A/) 2 f Af f Af \H(f - f')\ 2 dfdf> and hence that £ ~ 5f/Af. 



4.4 Multitaper Estimates 

While the lag window estimator is based on the idea that nearby frequencies 
provide independent estimates, the estimation is not very systematic, since, one 
should be able to explicity decorrelate nearby frequencies from the knowledge of 
the correlations introduced by a finite window size. This is achieved in multitaper 
spectral estimation. The basic idea of multitaper spectral estimation is to average 
the spectral estimates from several orthogonal tapers. The orthogonality of the 
tapers ensures that the estimates are uncorrelated for large samples (consider 
substituting h 1 (t)h 2 (t) for h{t) 2 in equation [17]). A critical question is the choice 
of a set of orthogonal tapers. A natural choice are the discrete prolate spheroidal 
sequences (dpss) or Slepian sequences, which are defined by the property that 
they are maximally localised in frequency. The dpss tapers maximize the spectral 
concentration defined as; 

\ = I W w\H(f)\ 2 df m] 
I-oo\H(f)\ 2 df 1 ] 
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Where in the time domain h(t) is strictly confined to the interval [0,T]. 

For given values of W and T there are a finite number of tapers which have 
concentrations (A) close to one, and therefore have well controlled bias. This 
number is known as the Shannon number and is 2WT. This sets an upper limit 
on the number of independent estimates that can be obtained for a given amount 
of spectral smoothing. 

A direct multitaper estimate of the spectrum is defined in equation |25|. 

/ MT (/)4?4 D (/) (25) 

The eigenspectra I® are direct spectral estimates based on tapering the data 
with the k th dpss function. As previously trial averaging can be included by using 
I DT rather than I D . More sophisticated estimates involve adaptive (rather than 
constant) weighting of the data tapers flPercival and Walden, 1993[ ). Multitaper 
spectral estimation has been recently shown to be useful for analysing neurobio- 
logical time series, both continuous processes ( Mitra and Pesaran, 1999|) and spike 
trains ( [Pesaran et al., 2000| ). 



4.4.1 Bias 



The bias for the multitaper estimate is given by equation 15 but with \H( 



replaced by an average over tapers 4 J2k=o 



4.4.2 Asymptotic Variance 

The asymptotic variance of the multitaper estimator, including trial averaging, is 
given by equation ^| . 



var{I^(f)} 



1 

N T K 



E{I (f)} 



(26) 



4.5 Degrees of freedom 

At this point it is useful to introduce the concept of the degrees of freedom (uq) 
of an estimate. The degrees of freedom is twice the number of independent esti- 
mates of the spectrum. Degrees of freedom is a useful concept as it permits the 
expressions for the variance of the different estimators to be written in a common 
format. 

var{I (/)} = (27) 

Where, 
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2N T K 



Degrees of freedom is also a useful framework in which to cast both finite size 
corrections and the confidence limits for the spectra and coherence. 

The variance of estimators of the spectrum can be estimated using internal 
methods such as the bootstrap or jackknife ( [Efron and Tibshirani, 1993| ), ( [Thom- 
son and Chave, 1991|) . Jackknife estimates can be constructed over trials or over 
tapers. If z/ is large (> 20), then the theoretical and Jackknife variance are in 
close agreement. If distributional assumptions can be validly made about the 
point process, theoretical error bars have an important advantage over internal 
estimates since they enable the understanding of different factors which enter into 
the variance in order to guide experimental design and data analysis. However 
Jackknife estimates are less sensitive to failures in distributional assumptions, and 
this provides them with statistical robustness. 

It is conventional to display spectra on a log scale. This is because taking 
the log of the spectrum stabilizes the variance and leads a distribution which is 
approximately Gaussian. 



4.6 Confidence intervals 



The expected values of the estimators and also their variance have been discussed 
for several different spectral estimators but it is desirable to put confidence inter- 
vals on the spectral estimates rather than standard deviations. 

As previously mentioned in section POI the averaging of direct spectral esti- 



mates from different trials yields, in the large sample limit, estimates which are 
distributed as x\n • m general for the other estimates a well known approxima- 



tion ( [Percival and Walden, 19^3 ) is to assume that the estimate is distributed as 



XZ ■ Confidence intervals can then by placed on estimates on the basis of this x 



I'D 



distribution. The confidence interval applies for the population spectrum S(f) 
and is obtained from the following argument. 



P 



qi < Xu < 92 



2p 



(2f 



Where P indicates probability, qi is such that P[xl < Qi] — P an d q% is such that 
P[xl > 92 ] = P- It follows that, 



qi<VoI X (f)/S(f)<q 2 =l-2p 



(29) 



Hence an approximate 100% x (1 — 2p) confidence interval for S(f) is given by, 



v I*(f)/q2<S(f)<v I*(f)/q 



(30) 
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For large uq (> 20) these confidence intervals do not differ substantially from 
those based on a Gaussian (±2 standard deviations) but at small u the difference 
can be substantial as for these values the xt n distribution has long tails. 



4.7 High Frequency limit 

The population spectrum goes to a constant value equal to the rate A in the high 
frequency limit. In practice spectra calculated from a finite sample will go to a 
value close to A but fluctuations in the number of spikes in the interval will lead 
to an error in this estimate. For a given sample the spectrum will go the value 



given by equation [31 . 



i K-l N T N n {T) 

nf - oo) = — e e e h k ( t jr (si) 

iV T^ fc=0n=l j=l 

Where f- is the j th spike in the n th trial and N n (T) is the total number of spikes 
in the n th trial. In the case of direct and lag window estimators the averaging 
over tapers need not be performed. 

This expression yields a value which is typically very close to the sample esti- 
mate of the mean rate[|. It is significant departures from this high frequency limit 
which are of interest when interpreting the spectrum as these indicate enhance- 
ment or suppression relative to a homogeneous Poisson process. 

4.8 Choice of estimator, taper and lag window 

The preceding section discussed the large sample statistical properties of direct, lag 
window and multitaper estimates of the spectrum. The choice of which estimator 
to use remains a contentious one ( Percival and Walden, 1993|) . The multitaper 



method is the most systematic of the estimators but the lag window estimators 
should perform almost as well for those spike train spectra which have reasonably 
small dynamic ranges[|. However, it is possible to construct spike trains with 
widely different time scales, which can possess a large dynamic range. In addition, 
the multitaper technique leads to a simple jackknife procedure by leaving out one 
data taper in turn. A further important property of the multitaper estimator is 
that it gives more weight to events at the edges of the time interval and thus 
ameliorates the arbitrary downweighting of the edges of the data introduced by 
single tapers. 

8 It is exactly the sample estimate of the mean rate for a rectangular taper. 
9 Dynamic range is a measure of the variation in the spectrum as a function of frequency and 
is defined as 10log w (^j§^). 



14 



If using the lag window estimator there are many choices available for both the 
taper and the lag window. The choice of taper is generally not critical provided 
that the taper goes smoothly to zero at the start and end of the interval. A 
rectangular taper has particularly large sidelobes in the frequency domain which 
can lead to significant bias. The choice of lag window is also usually not critical 
and typically a Gaussian kernel will be satisfactory. 



5 Estimating the Coherency 

Sample coherency, which may be estimated using equation |32], may be evaluated 
using any of the previously mentioned spectral estimators. The principle difference 
is that the direct estimator, in terms of which the other estimators are expressed, 
is given by equations |33] and [33] rather than |^ and [| 



C(f) 

JaU) 



J 12 



J?{f)J?(f) m 

h(t)e- 2 ™ ft dN a (t) 



(32) 
(33) 
(34) 



Where the N±(t) and A/" 2 (t) are simultaneously recorded spike trains from two 
different cells and X denotes the type of spectral estimator. Possible choices of 
estimator X include; D direct, DT trial averaged direct, LW lag window or MT 
multitaper. 



Lag window and multitaper coherency estimates may be constructed by sub- 
stituting I^(-) in place of I D (-) in equations [L8| and[25|. The estimates are biased 
over a frequency range equal to the width of the smoothing although the exact 
form for the bias is difficult to evaluate. 



5.1 Confidence limits for the Coherence 

The treatment of error bars is somewhat different between the spectrum and the 
coherency, since the coherency is a complex quantity. Usually one is interested is 
in establishing whether there is significant coherence in a given frequency band. In 
order to do this the sample coherence should be tested against the null hypothesis 
of zero population coherence. The distribution of the sample coherence under this 
null hypothesis is given below. 

P(\C\) = (v -2)\C\(l-\C\ 2 ) {uo/2 - 2) 0<|C|<1 (35) 

A derivation of this result is given in ( |Hannan, 1970| ). In outline the method is 
to rewrite the coherence in such a way that it is equivalent to a multiple correlation 
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coefficient flAnderson, 1984j ) . The distribution of a multiple correlation coeffient is 



then a known result from multivariate statistics. In the case of coherence estimates 
based on lag window estimators the appropriate vq may be used although this is 
only approximately valid because this method of derivation assumes integer vq/2. 
It is straightforward to calculate a confidence level based on this distribution. 



The coherence will only exceed yl — p l /( u o/2-i) i n p X 100% of experiments. In 
addition it is notable that the quantity (u Q /2 — 1)|C| 2 /(1 — \C\ 2 ) is distributed 
as ^2,1/0-2) under this null hypothesis. It is useful to apply a transformation to 
the coherence before plotting it which aids in the assessment of significance. The 



variable q = y — (z^o — 2)log(l — \C\ 2 ) has a Raleigh distribution which has density 

p(q) = qe~ q2 ^ 2 . This density function does not depend on vq and furthermore has a 
tail which closely resembles a Gaussian. For certain values of a fitting parameter^ 
j3, a further linear transformation r = /3(q— 0) leads to a distribution which closely 
resembles a standard normal Gaussian for r > 2. This means that for r > 2 one 
can interpret r as the number of standard deviations by which the coherence 
exceeds that expected under the null hypothesis. 

5.2 Confidence Limits for the Phase of the Coherency 

If there is no population coherency then the phase of the sample coherency is 
distributed uniformly. If, however, there is population coherency then the distri- 
bution of the sample phase is approximately Gaussian provided that the tails of 
the Gaussian do not extend beyond a width 2ti. An approximate 95% confidence 
interval for the phase QRosenberg et al., 1989| ; [Brillinger, 1974| ) is given below. 



(36) 



Where 4>(f), the sample estimate of the coherency phase, is evaluated using 
tan- 1 {/m(C)/ J Re(C)}. 



6 Finite Size Effects 

In the preceding sections error bars were given for estimators of the spectrum and 
the coherence. However these error bars were based on large sample sizes (they 
apply asymptotically as T — > oo). Neurophysiological data are not collected in 
this regime and, particularly in awake behaving studies where data is often sparse, 
corrections arising at small T are potentially important. In order to estimate 
the size of these corrections a particular model for the point process is required. 

10 A reasonable choice for j3 is 23/20. 
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The model studied was chosen primarily for its analytical tractability while still 
maintaining a non-trivial spectrum. 

The model and the final results will be presented here but the details of the 
analysis are reserved until appendix [A]. The model is a doubly stochastic inho- 
mogeneous Poisson process with a Gaussian rate function. A specific realization 
of a spike train is generated from the model in the following manner. Firstly a 
population spectrum Sc{f) is specified. From this a realization of a zero mean 
Gaussian process Xo{f) is generated. A constant A, the mean rate, is then added 
to this realization. This function is then considered to be the rate function for 
an inhomogeneous Poisson process. A realization of this inhomogeneous Poisson 
process is then generated. A schematic of the model is shown in figure |5|. 




Figure 5: Schematic illustrating the model for which finite size corrections to the 
asymptotic error bars will be evaluated, (a) A spectrum Sc(f) is defined (b) A 
realization \a(t) is drawn from a Gaussian process with this spectrum (c) The 
mean rate A is added to Ag?(£) to obtain \(t) (d) This rate function is used to 
generate a realization of an inhomogeneous Poisson process yielding a set of spike 
times. 

Technically this is not a valid process because the rate function A(£) may be 
negative. However, if the area underneath the spectrum is small enough then the 
fluctuations about the mean rate seldom cross zero and corrections due to this 
effect are negligible. In addition large violations of this area constraint have been 
tested by Monte Carlo simulation and the results still apply to a good approxi- 
mation. 

An important feature of this model is that the population spectrum of the spike 
trains is simply the spectrum of the inhomogeneous Poisson process rate function 
plus an offset equal to the mean rateQ The spectrum of the rate function is a 
positive real quantity and therefore within this model the population spectrum 
cannot be less than the mean rate at any frequency. Intuitively, the reason for this 
is that the process must be more variable than a homogeneous Poisson process at 
all frequencies. 

To make the nature of the result clear a simplified version is given in equation 
11 This result does not depend on the Gaussian assumption. 
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[37J This version is for the particular case of a homogeneous Poisson process (which 
has a flat population spectrum) and a rectangular tapeiQ- 



var{I x (f)} = A 2 
Where A is the mean rate. 



1 



(37) 



A sample based estimate of N?T\ is the total number of spikes over all trials. 
It is to be noted that finite size effects reduce the degrees of freedom. This result 
implies that there is a point beyond which additional smoothing does not decrease 
the variance further and this point is approximately when v is equal to twice the 
total number of spikes. The full result is given in equations |38] - |43|. 



nr{I*(f)} = 2E ™ 2 (38) 
«/(/) u ^ 2TN T E{I x {f)Y 1 ) 



Where, 



ft fitfdt If X = LW, D or DT 

C £ = \ j^j:jofk(t) 2 f k >(tydt If X = MT (40) 



k,k' 



f(t/T) = VTh(t) (41) 

and, 

$(/) = X hf + A[E{I X U)} - X hf ] + 2[E{I X (0)} - A*,] + [E{I x (2f)} - X hf ] (42) 

X hf = E{I x (f oo)} (43) 

C x is a constant of order unity which depends on the taper. When a taper is 
used to control bias some of the spikes are effectively disregarded and this has an 
effect on the size of the correction. The function f(t) has the same form as the 
taper h(t) but is defined for the interval [0, 1]. C x is the integral of the fourth 
power of / and obtains its minimum value of one for a rectangular taper. It is 
typically between 1 and 2 for other tapers. In the multitaper case cross terms 
between tapers are included. 



Equation |42] describes how the finite size correction depends on the structure 



of the spectrum. $(/) is the sum of four terms. The first term is the only 

12 The expression also holds approximately for the multitaper estimate provided all tapers up 
to the Shannon limit are used. 
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term which is present for a flat spectrum. The second term is a correction which 
depends on the spectrum at the frequency being considered. The next two terms 
depend on the spectrum at zero frequency and the spectrum at twice the frequency 
being considered. The latter three terms all depend on the difference between the 
spectrum at some frequency and the high frequency limit. Equation |42| applies 
provided that the spike train is well described by the model. However, this is not 
necessarily the case and a suppression of the spectrum, which cannot be described 
by the model, often occurs at low frequencies^. In the event that there is a 
significant suppression of the spectrum $(/) may become small or even negative. 
To avoid this a modified form for <&(/) may be used which prevents this. 

$(/) = A fe/ + 4max([E{I x (f)}-A hf ],0) + 2max([E{I x (0)} - A hf ], 0)... 

+ max([E{I x (2f)}-A hf ],0) (44) 

The above modification to the result is somewhat ad hoc so Monte Carlo sim- 
ulations of spike trains with enforced refractory periods have been performed to 
test its validity. These simulations demonstrated that, although the correction de- 
rived using H3 was significantly different from that obtained from the Monte Carlo 



simulations in the region of the suppression, equation 44 provided a pessimistic 



estimate in all cases studied. This increases confidence that applying finite size 



corrections using equation ^4| will provide reasonable error bars for small samples. 

Equation |39] gives the finite size correction in terms of a reduction in vq. The 
new u(f) may be used to put confidence intervals on the results, as described in 
section |4lj although the accuracy of the xl assumption will be reduced. In the 
case of the coherence an indication of the correction to the confidence level can 
be obtained by using the smaller of the two v(f) from the spike train spectra to 
calculate the confidence level using equation |35]. In all cases if the effect being 
observed only achieves significance by an amount which is of the same order as 
the finite size correction then it is recommended that more data be collected. 



7 Experimental Design 

Often it is useful to know in advance how many trials or how long a time interval 
one needs in order to resolve features of a certain size in the spectrum or the 
coherence. To do this one needs to estimate the asymptotic degrees of freedom 
uq. This depends on the size of feature to be resolved a, the significance level 
for which confidence intervals will be calculated p and the fraction of experiments 
which will achieve significance V. In addition the reduction in the degrees of 
freedom due to finite size effect depends on the total number of spikes N s and also 
Ch (see section |6]). 

13 Note that any spike train spectra displaying significant suppression below the mean firing 
rate can immediately rule out the inhomogeneous Poisson process model. 
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An estimate of vq may be obtained in two stages. Firstly a,p and V are 
specified and used to calculate a degrees of freedom v. Secondly the asymptotic 
degrees of freedom z/ is estimated using u, N s and Ch- The feature size a = 
(5 — A)/A is the minimum size of feature which the experimenter is content to 
resolve. For example, a value of 0.5 indicates that where the population spectrum 
exceeds 1.5 A the feature will be resolved. The significance level should be set to 
the same value that will be used for calculating the confidence interval for the 
spectrum, typically be 0.05. For a given p there is some probability V that an 
experiment will achieve significance. To calculate v one begins with a guess v g . 



Then qi is chosen such that P 



X 2 



> 



p/2. On the basis of this one then 



evaluates V, = 1 — $ [91/(1 + a)] where $ is the cumulative xl g distribution^. If 
V, is equal to the specified fraction V then v 
This procedure is readily implemented as a minimization of (V — V,{v g )) 2 on a 
computer. Having obtained v one can estimate v§ using the following expression. 



v g otherwise a different v g is chosen. 



1 



1 C h [l + 4a] 



v 



(45) 



2N S [1 + a\ 

Where the 4a is omitted from the numerator if a < 0. 

Figure ^ illustrates example design curves generated using this method. These 
curves show the asymptotic degrees of freedom as a function of feature size for 
different total numbers of spikes. 





Figure 6: Example design curves for the case when p = 0.05, V = 0.5 and 
Ch = 1.5. The three curves correspond to N s = 00 (solid), N s = 100 (dashed), 
N s = 50 (dotted). 



The existence of a region bounded by vertical asymptotes implies that as long 
as the total number of measured spikes is finite, modulations in the spectrum 



14 These formulae apply for a > 0. If a < then P 
should be used. 



< 



= p/2 and V, = $ + a)] 
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below a certain level cannot be detected no matter how much the spectrum is 
smoothed. These curves may be used to design experiments capable of resolving 
spectral features of a certain size. 

In the case of the coherence one calculates how many degrees of freedom are 
required for the confidence line to lie at a certain level as described in section |57L| . 



8 Line Spectra 

One of the assumptions underlying the estimation of spectra is that the population 
spectrum varies slowly over the smoothing width (W for multitaper estimators). 
While this is often the case there are situations in which the spectrum contains 
very sharp features which are better approximated by lines than by a continuous 
spectrum. This corresponds to periodic modulations of the underlying rate, such 
as when a periodic stimulus train is presented. In such situations it is useful to 
be able to test for the presence of a line in a background of colored noise (i.e. in 
a locally smooth but otherwise arbitrary continuous population spectrum). Such 
a test has been previously developed, in the context of multitaper estimation, for 
continuous processes (Thomson, 1982[ ) and in the following section the analogous 
development for point processes is presented. 



8.1 F-test for point processes 

A line in the spectrum has an exactly defined frequency and consequently the 
process N(t) has a non-zero first moment. The natural model in the case of a 



single line is given by equation 46 



E{dN(t)}/dt = A + Aicos(2 7 r/ 1 t + 0) (46) 

A zero mean process (N) may be constructed by subtraction of an estimate 
of Xot. Provided that the product of the line frequency (/i) and the sample dura- 
tion^) is much greater than one the sample quantity N(T) /T is an approximately 
unbiased estimate of A . The resultant zero mean process N has a Fourier trans- 
form which has a non-zero expectation. 



Jk(f) = / h k (t)e- 2mft dN(t) (47) 

J — oo 

E{J k (f)} = c x H k {f-h)+<*H k {f + f x ) (48) 

Where, 

Cx = Aie^/2 (49) 
In the case where / > and fi > W, 

E{J k (f)} ~ Cl H k {f - fx) (50) 
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The estimates of Jk(fi) from different tapers provide a set of uncorrelated 
estimates of CiH k (0). It is hence possible to estimate the value of C\ by complex 
regression. 

Cl = Ejmm (51) 

Under the null hypothesis that there is no line in the spectrum (c± = 0) it may 
readily be shown that E{&i} = and var{di} = S(fi)/ J2k l-^fc(0)| 2 . The residual 
spectrumQ, which has the line removed, may be estimated using equation [52. 



$(/) = -?£l-W)- W/-/0I 3 (52) 



K k 



In the large sample limit the distributions of both c% and S(f\) are known 
( Percival and Walden, 1993|) and may be used to derive relation |k| 



^2,2{K-l) (53) 



Where = denotes 'is distributed as'. 

The null hypothesis may be tested using this relation and, if rejected, the 



line can be removed using equation [52| to estimate the residual spectrum. It is 
worth noting that although relation ^ was derived for large samples the test is 
remarkably robust as the sample size is decreased. Numerical tests indicate that 
the tail of the F distribution is well reproduced even in situations where there are 
as low as 5 spikes in total. 



8.2 Periodic Stimulation 

A common paradigm in neurobiology where line spectra are particularly important 
is that of periodic stimulation. When a neuron is driven by a periodic stimulation 
of frequency fi the spectrum may contain lines at any of the harmonics nf\. 
Provided that j\ > 2W the analysis of section |8.1| applies with each harmonic 
being separately tested for significance. 

The first moment of the process, which has period 1/fi, is given by equation 
|5"3| and may be estimated using c n . 

K*) = A o + E A n cos(2vrn/ 1 t + <f> n ) (54) 

n 

Where A n = 2|c n |, (fi n = tan _1 {/m(c n )/i?e(c n )}, the sum is taken over all the 
significant coefficients. 

15 It is also possible to estimate a residual coherency. In order to do this one uses a residual 
cross-spectrum S xy (f) = ± Y,k( J k(f) ~ ^iMf ~ /i))*(W) " c\H k {f - fx)), together with 
the residual spectra to evaluate the usual expression for coherency. 
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This rate function X(t) is the average response to a single stimulus or impulse 
response. The coefficients the Fourier series representation of X(t). 



8.3 Error Bars 

It is possible to put confidence intervals on both the modulus and the phase 
of the coefficients c n . For large samples(> 10 spikes) the real and imaginary 
parts of c n are distributed as independent Gaussians each with standard deviation 

!. For c n /a n > 3 the distribution of \c n \ is well 



a n = ^S{nf x )/{2T, k \H hK 
approximated by a Gaussian centered on \c n \ and with standard deviation o n . In 
addition the estimated phase angle (4> n ) is also almost Gaussian with mean <p n and 
standard deviation cr n /\c n \. Approximate error bars or confidence intervals may 

be obtained using a sample based estimate of a n , cr n = \J S{nfi) / {2 J2k \Hk(0)\ 2 )- 
Estimating error bars for the impulse response function is more involved due 
to their non-local nature (if one of the Fourier coefficients is varied the impulse 
response function changes everywhere). It is therefore of interest to estimate a 
global confidence interval, defined as any interval such that the probability of the 
function crossing the interval anywhere is some predefined probability. A method 
for estimating a global confidence band is detailed in (|Sun and Loader, 1994j) and 
outlined here. First a basis vector $(£) is constructed. 



&icos(2irfit) 

a N cos(2%f N t) 
aisin(27r fit) 



(55) 



a n sin{2n f ^t) 

Where N is the total number of harmonics. 

The elements of this vector have unit variance and a standard approximation 
may be applied. 

P{sup\X{t) - E{\{t)}\ > c||$(t)||) < 2(1 - N{c)) + {k/ii)e~ c2 ' 2 (56) 

Where sup is the maximum value of its operand, | |$(t) 1 1 denotes the length of vec- 
tor $(£), N(c) is the cumulative standard normal distribution and k is a constant. 
k may be evaluated by constructing the 2 x N matrix X(t) = d<f>(t) j dt], 
forming its QR decomposition ( [Press et al., 19TJ2" ) and then evaluating k = 



J, 1 \R 22 (t)/R n (t)\dt. 

Confidence intervals for the residual spectrum are calculated in the usual man- 
ner (using xt) although at the line frequencies the interval is slightly broadened 



due to the loss of a degree of freedom incurred by estimation of c n . Section |H 
contains an example application of the methods described in this section. 
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9 Example Spectra 



Figure |7| is a spectrum calculated from data collected from a single cell recorded 
from area PRR in the parietal cortex of an awake behaving monkey during a de- 
layed memory reach task ( [Snyder et al., 1997]) . The spectrum is calculated over 
an interval of 0.5s during which the firing rate is reasonably stationary and is av- 
eraged over 5 trials. The spectrum shows two features which achieve significance. 
There is enhancement of the spectrum in the frequency band 20-40 Hz indicating 
the presence of an underlying broad band oscillatory mode in the neuronal firing 
rate. In addition there is suppression of the spectrum at low frequencies. As dis- 
cussed previously a suppression of this sort is consistent with an effective refactory 
period during which the neuron is less likely to fire. Care must be taken at low 
frequencies since at frequencies comparable to the smoothing width the spectrum 
is particularly sensitive to any non-stationarity in the data. 




Figure 7: (a) Gaussian kernel (100 ms width) smoothed firing rate with la error 
bars based on a stationarity assumption. The vertical lines indicate the period 
over which the spectrum was calculated. A light is flashed at time zero and the 
spectrum is evaluated over the interval when the monkey is required to remember 
the target location, (b) The spectrum evaluated over this interval using a lag 
window estimator with a 40% cosine taper and a Gaussian lag window of width 
3.5 Hz. 95% confidence limits are shown with the finite size correction included 
(this typically resulted in a decrease in v(f) from about 50 to 36). The horizontal 
line indicates the high frequency limit. (c) The same spectrum evaluated using a 
multitaper estimator. A bandwidth (W) of 5 Hz was used allowing 5 tapers. Both 
estimators have the same degrees of freedom. 
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10 Example Coherency 



To illustrate the estimation of coherency simulated spike trains were generated 
from a coupled doubly stochastic Poisson process. For a given trial a pair of rate 
functions were drawn from a Gaussian process. The realizations share a coherent 
mode which is linearly mixed into the rates of both cells. These coupled rate 
functions are then used to independently draw a realization of an inhomogeneous 
Poisson process for each cell. Using this method 15 trials of duration 0.5s were 
generated. The coherent mode was set such that the population coherence was 
a Gaussian of height 0.35 and standard deviation 5 Hz centered on 20 Hz. The 
phase of this mode was set to 180°. Figure ^indicates that this coherent mode is 
reasonably estimated. 



1 4 




Figure 8: (a) Coherence (left axis) and the phase of the coherency (right axis). 
Fifteen trials of 0.5s duration were simulated using a doubly stochastic Poisson 
process as described in the text. A multitaper estimator with a smoothing width of 
7 Hz was used. Finite size corrections were used and resulted in 25% reduction in 
the degrees of freedom. A horizontal line has been drawn at the 95% confidence 
level under the null hypothesis of no coherency. Where the null hypothesis is 
rejected the phase of the coherency is estimated and shown with an approximate 
95% confidence interval, (b) The standardized coherence is a transformation which 
maps the null distribution onto an approximately standard normal variate (as 
described in section |5.1| ). The estimated coherence at 20 Hz would therefore lie 
at three standard deviations if there were no population coherence. 



11 Example Periodic Stimulation 

An example of an analysis of a periodic stimulus paradigm is shown in figure || 
The data is a single cell recording collected from the barrel cortex of an awake 
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behaving rat during periodic whisker stimulation at 5.5 Hz ( |Sachdev et al., 1999| ). 
There is a single trial of duration 50s. 




Figure 9: Response to a periodic stimulation of frequency 5.5 Hz. (a) Impulse 
response function with global 95% confidence interval (b) \c n \ versus index n with 
95% confidence interval. Dots indicate points which achieved significance in the 
F-test. (c) Residual spectrum with finite size corrected confidence interval. A 
multitaper spectrum with 100 tapers and a bandwidth of 1.5 Hz was used initially 
to avoid overlap of harmonics. This spectrum was then further smoothed using a 
Gaussian lag window with standard deviation 9 Hz. (d) The coefficient phases 4> n 
(in radians) versus index n after subtraction of a fitted straight line of gradient 
27r/3 ± 0.01. The black dashed lines are a 95% confidence interval about zero. 

The estimated impulse response function \(t) is seen to have two distinct sharp 
peaks outside of which the response does not differ significantly from zero. The 
moduli of the Fourier coefficients are significant out ton = 25. This automatically 
sets the smoothing of \(t) as structure on a time scale of less than 1/(25 x 5.5) = 
7 ms does not achieve significance. Note that the coefficients are enhanced at 
multiples of 6 (i.e.~ 33 Hz) which comes from having two peaks in the time 
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J k a n U) = So h k (t)e-^dN2t) 

W) = J k a n u)J k ru) 



Table 1: The basic direct spectral estimator in terms of which the other estimators 
can be written. For clarity the superscript D on the direct spectral estimate has 
been omitted. The index n labels trials, index k labels tapers, and indices a and 
b label cells. 

domain X(t) which are separated by ~ 30 ms. The phase of the coefficients closely 
follows a straight line but there is a small periodic deviation from this line which is 
again at index multiples of 6. The gradient of the straight line depends on the time 
delay of the response. The residual spectrum was calculated by first evaluating 
a multitaper estimate from which the significant harmonics were removed. This 
spectrum had a bandwidth of 1.5 Hz chosen to avoid overlap of the harmonics 
leading to the multitaper estimate being undersmoothed. A further smoothing 
was performed using a lag windowQ and the resultant spectrum, displays a slight 
but significant suppression relative to a Poisson process out to almost 200 Hz. 
Such a spectrum is characteristic of a short time scale refractive period. The 
residual spectrum is particularly useful because rate non-stationarity has been 
removed. 

12 Summary 

It is the belief of the authors that spectral analysis is a fruitful and under exploited 
analysis technique for spike trains. In this paper an attempt has been made to 
collect the machinery necessary for performing spectral analysis on spike train 
data into a single document. Starting from the population definitions the statis- 
tical properties of estimators of the spectrum and coherency have been reviewed. 
Estimation methods for both continuous spectra and spectra which contain lines 
have been included. In addition new corrections to asymptotic error bars have 
been presented which increase confidence in applying spectral techniques in prac- 
tical situations where data is often sparse. Tables 1 to 5 summarize the important 
formulae. Matlab software implementing the methods discussed in this paper is 
available from xxx.lanl.gov/archive/neuro-sys. 

16 The previous theory developed for lag window estimators applies to this hybrid esti- 
mator with |-ff(-)| 2 replaced by J2k=o Wk(;)\ 2 in equation ^ and |H(-)| 2 replaced by 
Ti Hk=o l^fe(-)! 2 in equation ||. 
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Table 2: The different estimators and the large sample degrees of freedom uq of 
estimates of the spectrum (ab = 11). The indices on the 1^ are as follows, ab 
label the cells from which the estimates are constructed. The index k labels the 
taper and n labels the trial. 





Equation Eq. 


Comment 


Variance 




use fo for asymptotic 
or if using finite 
size correction 


Degrees 
of freedom 


1 _ 1 1 57J 


See text for definitions 
ofCf and $(/) 


Confidence 
(1 - 2p) x 100% 


[^/ x (/)/? 2 ,^ x (/)/gi] a 


9i s.t P[x^ < gi] = p 
q 2 s.t P[x^ > g 2 ] = p 



Table 3: Main formulae required for estimating spectral error bars. Refer to 
section § for additional information. 





Equation Eq. 


Comment 


Coherency 


C X (f) = ~fa 

V aa bb 


32 


1 




Distribution 
for coherence 


P(\C\) = — 2)|C|(1 — |C7| 2 )t-/2-2) 


35 


1 


Under null 
hypothesis 7 = 


Confidence 
for phase 




36 


1 


Approx. 

95% 



Table 4: Main formulae required for coherency estimation. Refer to section |5| for 
additional information. 
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Equation Eq. 


Comment 


Complex amplitude 
of line 


Cl £ fc |**(0)|» 


51] 




F-test to access the 
significance of a line 


^JJ fc (/i) eiH fe (o)P H2(K-1) I 


53] 


Null 

ci = 


Residual spectrum 




52] 





Table 5: Main formulae required for the detection and removal of a line from the 
spectrum. Refer to section |8] for additional information. 

A Derivation of Finite Size Correction 

The following is an outline derivation of the finite size corrections described in 
section ||. Firstly the characteristic functionals flBartlett, 19661 ) f° r the processes 
N and the inhomogeneous Poisson process rate X(t) are related. 

Cjf(6(t)) = E{exp(i F 6{t)dN)} = E x {exp( f \(t)b(6(t))dt)} (57) 
Jo Jo 

b(G(t)) = exp[i0(t) - - C 6(t')dt'] (58) 

T jo 

Under the Gaussian process assumption for \(t) this integral may be done. 

C w = exp[- ( ( b(t)k(t,t')b(t')dtdt' + A / b(t)dt] (59) 
2 jo J o Jo 

A{t,t') = E x {(X(t) - A)(A(t') - A)} (60) 

Note that A denotes the mean rate. Taking the log of the characteristic functionals 
now yields the following relation between the resultant cumulant functionals. 

K w =lnE{exp(i [ 6(t)dN)} = - [ [ b(t)A(t,t')b(t')dtdt' + A / b(t)dt (61) 
Jo 2 Jo Jo Jo 

Next 6 {t) is chosen appropriately and substituted into Kjj. The form for 9(t) 
which is required to obtain the covariance of multitaper estimators is; 



id(t) = ^i/i fc (t)e~ 2m/l * + e 2 h k (t)e 2mflt + 6 3 h k ,(t)e- 2mf2t + 6 4 h k ,(t)e 27Tif2t (62) 
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Substituting into the cumulant functional for N yields; 



Kjf = lnE{exp(6 1 J^(f 1 ) + fcJ?*(/i) + 9 3 J°(f 2 ) + 9 4 J°*(f 2 ))} (63) 

Where J® is the Fourier transform of the data tapered by a function indexed by 
k. Application of the cumulant expansion theorem ( [Ma, 198"5|) then leads to; 



K w = E{exp(9 1 J^(f 1 ) + MH/i) + 9 * J k>{h) + Mj£*(/ 8 )) ~ (64) 
This may then be differentiated and set to zero. 



K, 



Imno 



d6[defdB%dei 



E{Jk\h)Jk m *{h)J§ n {h)J^r{h)} c 



jDm* 



Dn, 



tDo* I 



01=02=8 3 -- 



(65) 

Moments of the estimators may be expressed in terms of these cumulant deriva- 
tives. The expressions are simplified by the fact that all cumulant derivatives 
which have indices summing to an odd number are zero because N is a zero mean 
process. 

E{I°(f)} = K 1W0 (66) 



var{I MT (f)} 4F^«/)»^(/)} 

" k=0 fc'=0 



cov{I°(f), /£(/)} = K W10 K 0W1 + Kmi + K W01 K 01 



rD 



10 



(67) 
(68) 



The problem has now been reduced to that of calculating these derivatives within 
the model. This is done by substituting the expression for 8(t) into the RHS of 
equation Considerable algebra then leads to the following exact result. 



K, - K A 4- K 

J-^-lmno — Iy lmno ' IV I 



B 

Imno 



(69) 



Where, 



1 ^ l\m\n\o\ 
/ ""°^2 ; ^ lI/,!lI///,!lI,/,!lIo,! 

li,rm,Tii,oi % i % i 



-Hxjh 

T 



■Hi(f: 







\-Hi(f 2 y] 


T 




T 



T 

02+04 



m2+ni4 



1 l 3 ,m 3 ,n 3 ,o 3 



(70) 



Where J2i h — I an d cases where li + 1 2 — I or l 3 + 1± = I are excluded (and also 
for n, m, o). 

ttSS = / S x (f)H h+mi+ni+0l {Mh - rm) + f 2 (m - 0l ) - /] 

^; +m3+n3+03 [A(/3 - m 3 ) + f 2 (n 3 - o 3 ) - f]df (71) 
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Where S\(f) is the spectrum of the Gaussian process and Hi is; 

/oo 
h l (t)exp(-2mft)dt (72) 
-oo 

H (f) = Texp(-iTrfT)sinc(7rfT) (73) 

l m n o 

C = ^EEEDH? l ][' , ][8 ]iW- +s [/x(p - g) + h{r - s)\ ... 

p=0 g=0 r=0 s=0 



r-^i(/i) _ 


(i-p) 




(m-q) 


\-Hi{f»Y 


(n—r) 


\-Hi(h)*] 


T 




T 




T 




T 



(o-s) 

(74) 



The preceding result is somewhat cumbersome but readily evaluated computation- 
ally for a given spectrum. The expression simplifies greatly when only frequencies 
above the smoothing width are considered and many of the terms may be ne- 
glected. Restricting attention to the second order properties there are only a few 
remaining dominant terms. Terms from Kiooi lead to the previously discussed 
asymptotic results but there are corrections which arise from the term Kim- As- 
suming that the population spectrum varies slowly over the width of the tapers 
leads to the result given by equations ^ - |43|. The validity of this assumption has 
been tested computationally and was found to be very accurate even for spectra 
with sharp peaks. 
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